arXiv:1502.01046vl [physics.chem-ph] 3 Feb 2015 


Assessing the Performance of the Diffusion Monte Carlo Method as Applied to the 

Water Monomer, Dimer, and Hexamer. 

Joel D. Mallory]^ Sandra E. Brownjf] and Vladimir A. Mandelshtanjl] 

Department of Chemistry, University of California, Irvine, 

1102 Natural Sciences II Irvine, California 92691, USA 
(Dated: February 5, 2015) 

The Diffusion Monte Carlo (DMC) method is applied to the water monomer, dimer, and hexamer, 
using q-TIP4P/F, one of the most simple, empirical water models with flexible monomers. The bias 
in the time step (At) and population size (N w ) is investigated. For the binding energies, the bias in 
At cancels nearly completely, while a noticeable bias in N w still remains. However, for the isotope 
shift, (e.g, in the dimer binding energies between (^ 0)2 and ( 020 ) 2 ) the systematic errors in N w do 
cancel. Consequently, very accurate results for the latter (within ~ 0.01 kcal/mol) are obtained with 
relatively moderate numerical effort (N w ~ 10 3 ). For the water hexamer and its (D20)6 isotopomer 
the DMC results as a function of N w are examined for the cage and prism isomers. For a given 
isomer, the issue of the walker population leaking out of the corresponding basin of attraction is 
addressed by using appropriate geometric constraints. The population size bias for the hexamer is 
more severe, and in order to maintain accuracy similar to that of the dimer, the population size 
N w must be increased by about two orders of magnitude. Fortunately, when the energy difference 
between cage and prism is taken, the biases cancel, thereby reducing the systematic errors to within 
~ 0.01 kcal/mol when using a population of N w = 4.8 x 10° walkers. Consequently, a very accurate 
result for the isotope shift is also obtained. Notably, both the quantum and the isotope effects for 
the prism-cage energy difference are small. 


Introduction 

Diffusion Monte Carlo (DMC) jlMj has engendered sig¬ 
nificant attention in the literature because it is one of the 
few numerical methods that enable the computation of 
the ground state of many-body systems. Large and cum¬ 
bersome basis sets that scale exponentially with the sys¬ 
tem size and constitute an essential component of most 
variational methods are nonexistent in DMC[5j. Conse¬ 
quently, the computational cost of DMC increases rather 
slowly with the particle numberJB], which explains the 
method’s attractiveness for treatments of systems with 
many degrees of freedom. 

DMC is appealing for theoretical studies of condensed 
matter systems and has developed a reputation for its 
widespread application to weakly-bound ensembles of 
molecules that interact through hydrogen bonding and 
dispersion forces 0 [5]. Some examples of systems that 
have been examined with DMC are Bose condensates of 
parahydrogen and helium [9lU2j: sheets of graphite and 
diamond[T51; fee crystallized xenon at 0 K0; HF, HCN, 
and SF 6 trapped in clusters of argon and helium [51 ID]: 
and water clusters with a special emphasis placed on the 
water hexamer [T5f - f22| . 

As originally proposed by Anderson[l, 2] DMC takes 
advantage of the similarity between the imaginary time- 
dependent Schrodinger equation and the diffusion equa¬ 
tion. The method employs a population of random walk- 
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ers (or J-functions) that evolves in imaginary time r and 
samples the configuration space to collectively represent 
the lowest energy eigenstate of the system defined by a 
potential energy surface (PES). 

A longstanding issue with DMC is that it suffers from 
inherent sources of systematic errors that arise due to 
the use of finite values for the time step At and pop¬ 
ulation size N w , as well as the need to introduce a 
population control mechanism in the form of branching 
steps [71151 H01125] . We have found that an analysis of the 
behavior and extent of the biases is often absent or ig¬ 
nored altogether when the results from DMC simulations 
are reported (see, e.g., refs, fill H2l H5HT71 l20[l . Ideally, 
the information on the DMC energy estimate at finite 
values of At and N w can be analyzed and extrapolated 
to At —► 0 and N w —► 00 . Yet, ref. [9] suggests that the 
bias of the DMC energy estimate, specifically, the pop¬ 
ulation size bias, can be very nontrivial and as such is 
often difficult to extrapolate. 

Nevertheless, efforts have been underway for some time 
to substantially reduce or even eliminate the biases with 
respect to the time step and population size. The price 
to be paid, however, is that the algorithms are more com¬ 
plicated than the simple version developed by Anderson. 
These variants of DMC usually implement importance 
sampling (IS-DMC) whereby a drift velocity term is in¬ 
corporated into the imaginary time Schrodinger equa¬ 
tion to drive the random walkers into regions of config¬ 
uration space where the wavefunction is largest. This 
procedure often utilizes an optimized guiding or trial 
wavefunction [TSJ I25H2H] . 

The time step error was studied extensively by Umri- 
gar in ref. [25] It was determined that the bias could be 
considerably diminished if diffusive displacements of the 
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random walkers are accepted or rejected based on a ra¬ 
tio of values for the trial wavefunction at the current and 
previous points in configuration space. This methodology 
can only be implemented in conjunction with IS-DMC, 
as the probability for acceptance of the diffusion move is 
directly linked to a Gaussian distribution containing the 
aforementioned drift term. It was demonstrated that this 
algorithm does, in principle, reduce the time step error, 
but only if an accurate form for the trial wavefunction is 
known a vriori [25]. Thus, the most practical approach to 
remove the time step error is still to run several simula¬ 
tions at different values of At and extrapolate the result 
to At = 0. 

An “improved DMC method” proposed in ref. [23] es¬ 
tablishes a fixed population of random walkers with un¬ 
equal weights. Spurious correlations among the other¬ 
wise independent motions of the walkers contribute to 
the population size (or population control) bias and arise 
when the weights are reconfigured during the branching 
steps. Therefore, stochastic reconfigurations are carried 
out with the intention of minimizing the fluctuations in 
the weights as much as possible. This technique may 
decrease the population size bias and the statistical fluc¬ 
tuations in the energy [13i [23]. 

Another DMC variant, called norm-conserving DMC, 
exactly balances the flux of walkers entering and exiting 
the ensemble such that the population remains constant 
for all iterations of the simulation[7:. Weights for the 
walkers are not introduced except through a mean-field 
scheme that is invoked to decrease the population size 
bias, which seems to scale as ~ 1 /N w . It was asserted 
that this approach could entirely eliminate the popula¬ 
tion size bias, but in reality the bias will vanish strictly 
in the limit of N w —> oo HUM! [HI [27]. 

The issue of a non-vanishing population size bias was 
reported by Boninsegni[2] upon observing discrepancies 
between DMC energies for helium clusters and energies 
obtained through an alternative method [10|. In a subse¬ 
quent study of parahydrogen clusters, the bias was thor¬ 
oughly investigated and found to never completely dis¬ 
appear, even in an attempt to extrapolate the result to 
N w —>• oo. Indeed, the dependence of the energy on the 
reciprocal walker number proved to be nontrivial, and the 
resultant curve seemed to diverge in the limit of large 
walker numbers. In light of this study we believe that 
the population size bias must be thoroughly examined 
whenever DMC is utilized in treatments of many-body 
systems. The suggested divergence of the DMC energy 
at large N w should make all reported DMC results ques¬ 
tionable, whenever such a study is lacking. 

Furthermore, ref. [3| utilizes the DMC variant outlined 
in ref. [25] which requires IS-DMC such that the popula¬ 
tion size bias happens to be quite sensitive to the choice 
of the trial wavefunction. If a substantial discrepancy ex¬ 
ists between the trial and exact ground state wavefunc- 
tions, the systematic bias with respect to N w becomes ex¬ 
tremely pronounced. Since such a discrepancy is unavoid¬ 
able, the pathological convergence with respect to the 


number of random walkers may be an intrinsic property 
of IS-DMC. In other words, IS-DMC may not be ideal for 
treating complex systems with complicated ground state 
wavefunctions, as any attempt to optimize the trial wave- 
function by finding a suitable parametrizaion is bound to 
fail. Water clusters seem to provide a good example of 
such systems, in which the complexity is due to the ex¬ 
istence of two different time scales corresponding to the 
slow inter- and fast intra-molecular degrees of freedom. 

The water hexamer has been intensely studied from 
many different experimental and theoretical angles be¬ 
cause it represents the smallest cluster of water molecules 
that still maintains a three-dimensional configuration [3, 
EE5J EE3 l2Qj [22] . This fundamental building block of water 
and ice retains many of the physical properties of bulk 
water and the ability to exist in discrete isomeric forms. 
In order to avoid dealing with the two time scales most of 
the DMC simulations of water clusters and, in particular, 
the water hexamer have invoked the frozen monomer ap¬ 
proximation such that the intramolecular degrees of free¬ 
dom for the constituent water molecules were collectively 
neglected [I3HI31 [2D| [2Tj . This approach assumes that the 
intra- and inter-molecular motions can be adiabatically 
decoupled, and as such, lead to a substantial simplifica¬ 
tion of the water dynamics and its numerical treatment. 
Unfortunately, this approximation fails to describe many 
important properties of water, including the notable iso¬ 
tope effect due to the substitution of hydrogen atoms by 
deuterium atoms. 

Consequently, in the present study we consider a flex¬ 
ible water model in conjunction with full-dimensional 
DMC. Nonetheless, our focus is not to carry out accu¬ 
rate first-principles simulations of water that would, in 
particular, employ an “accurate”, possibly ab initio, wa¬ 
ter PES. Instead, we focus our attention on the method¬ 
ological issues in an attempt to assess the performance 
and various other aspects of DMC implementation for a 
complex many-body system—a prime example of which 
is the water hexamer. In order to make our goal feasible, 
we minimize the computational cost by employing the q- 
TIP4P/F PES; one of the least expensive empirical, flex¬ 
ible water models available. This potential has proved 
to be reasonable in describing thermodynamic properties 
of water with quantum nuclear effects included [28 i. We 
note that presumably more accurate PES’s exist, such as 
the WHBB[2§] or HBB2-pol[3^, and the most recent and 
most accurate MB-pol PES[3l]. which are all ab initio- 
based. Unfortunately, all these potentials are very ex¬ 
pensive, and by orders of magnitude more expensive than 
q-TIP4P/F. We believe that the present study using the 
inexpensive PES will allow us to carry out DMC simu¬ 
lations employing sufficiently long projection times and 
large population sizes to be able to assess the method¬ 
ology before it is applied to more realistic, albeit more 
expensive, PES’s. 

Notably, Bowman and co-workers [15] have already re¬ 
ported DMC results on the ground state energies of three 
isomers of the WHBB water hexamer (cage, prism, and 
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book). The authors briefly acknowledge in the supple¬ 
mental information that a strong dependence of the iso¬ 
mer energies on the population size is evident, but the 
extent and behavior of the bias was not rigorously char¬ 
acterized. When taking the energy difference, the large 
systematic errors may or may not cancel, and a nontriv¬ 
ial bias can arise for the energy difference with respect 
to the population size and/or time step. Most impor¬ 
tantly, however, the issue arising due to leakage of the 
random walkers out of their original basin of attraction 
has not been thoroughly addressed. We have observed 
that the previous DMC treatments of the water hexamer 
using rigid monomers [T71 [TH] 2S EL] either did not con¬ 
sider this problem, or did not consider it adequately, but 
we are unaware of whether it poses a serious problem 
for clusters comprised of rigid water monomers. For ex¬ 
ample, ref. l20l does mention this problem, but makes an 
impression that it is not serious. The authors of the latter 
work suggested that using a simple constraint based on 
monitoring only the energies of random walkers solves the 
problem. In the case of a flexible water hexamer (at least 
with the q-TIP4P/F PES), the migration of the popula¬ 
tion out of the basin of attraction in question occurs for 
all the isomers considered in our test calculations (some 
of which are not reported here). Therefore, a numeri¬ 
cally inexpensive and effective solution to this problem 
appears to be nontrivial. In particular, a simple energy- 
based constraint would fail, because the energy differ¬ 
ences between different isomers (and the corresponding 
quantum energies) are small. In contrast, the authors of 
ref. Q31 proposed the use of a relatively simple geometric 
constraint that seemed to prevent the population migra¬ 
tion in their DMC simulation of Ar„HF van der Waals 
clusters. Our approach involves implementation of more 
than one simple geometric constraint. 

In what follows, we present a comprehensive assess¬ 
ment of the performance of DMC as prescribed by 
Anderson [T] E and its sources of systematic bias in ap¬ 
plications to the water monomer, dimer, and hexamer. 
In addition, we consider the isotope substitution of hy¬ 
drogens by deuteriums and, consequently, investigate the 
isotope shifts and their dependence on both the popula¬ 
tion size and the time step. 

The Diffusion Monte Carlo Method 

Consider an 7V-particle system described by the Hamil¬ 
tonian 


fusion equation [HE 

= _(ff _ E ief )V = -Ttf - (V - £ re f)d> (2) 

OT 

with the constant shift E re f defining the energy reference. 
This equation can be solved and recast in terms of the 
imaginary time propagator 

W(r;r) = e (B - f " A)T ^(r, 0) (3) 

By expanding d' (0) in the eigenbasis of H and substitut¬ 
ing into Eq. 0, 

T(r;r) = ^ c fc e (gref ~ gfc)r y> fc (r) 
k 

= e^- E ^ T J2 c ke {Eo ~ Ek)T Mr) ( 4 ) 

k 

one can see that in the r —> oo limit T(r; r) becomes pro¬ 
portional to the ground state wavefunction ’Fo, as all the 
other contributions tend to zero exponentially, relative 
to the contribution of the ground state. Furthermore, 
Eq. @ also implies that 

f o, if E le { < E 0 

lim IjTlI = < oo, if E lei > E 0 (5) 

[ const, if E le f = E 0 

That is, the solution of Eq. ([0) is stationary only when 
E re { coincides with the ground state energy Eq. 

Eq. |2| can be solved numerically by introducing a fi¬ 
nite time step At and invoking the split-operator approx¬ 
imation for the short-time propagator 

T(r; t + At) » e {E ^~ V)AT e~ fAT ^{v- t ) (6) 

Assuming the wavefunction d' is non-negative at all po¬ 
sitions, it can be represented by an ensemble of delta 
functions WjS (r — r^) or random walkers = r^\r) 
(j = 1, ...,N W ) that evolve in time t together with their 
weights Wj = Wj{r). 

Operating with the kinetic energy propagator on the 
i-th component of the j -th d-function, 

3 N 

<5(r — r (j) ) = S (rj - r| j) ) 

»=l 

gives a Gaussian distribution that governs the diffusive 
displacements for a single component of the j-th random 
walker in configuration space 


3N „ 2 

H = T + V(r) = J2^- + V(r) ( 1 ) 

where V(r) defines the PES, m, are particle masses, and 
r = (ri, ...,r 3 jv) € R 3Ar is the coordinate vector. 

The DMC method of Anderson exploits the ap¬ 
parent isomorphism between the imaginary time- 
dependent Schrodinger equation and the standard dif- 


m i (r l - r ^) 2 

2Kt 

(7) 

Consequently, the kinetic energy part of the random 
walker propagation can be represented by shifting the 
i -th component of j-th walker randomly, 

r [ j \r + At) = r[ j \r) + z» 


exp 


A rpi 2 
2rrii 




( 8 ) 












4 


where z* is a Gaussian random number with standard 
deviation 



The action of the potential energy propagator on the 
walker population can be implemented in many different 
ways. The most straight-forward recipe is to multiply 
each weight Wj by the factor 


Pj = exp 


(£ ref -V(r«))Ar 


(9) 


However, this will soon make the random walker repre¬ 
sentation of the wavefunction inefficient due to the ap¬ 
pearance of random walkers with either very small or 
very large weights. Here we adapt the simplest version of 
DMC, in which the weights are all the same, Wj = 1, but 
upon each advancement by At a “branching” procedure 
is implemented, in which some walkers are replicated and 
some are killed. 

The value of of Pj computed upon each diffusion dis¬ 
placement determines the number of copies (n ? ) of the 
random walker (i.e. r^) that will be retained in the 
ensemble. Precisely, 


[Pj] + !» 
\Pj\ > 


if {Pj} > £ 
otherwise 


( 10 ) 


where £ € [0,1] is a uniformly distributed random num¬ 
ber. In particular, the random walker is killed when 
nj = 0 . 

Upon completion of each branching step the value E re { 
is updated to reflect changes in the population size N w 
according to 


where 


-Uref(T) 


V(t) 


N w (t)-N w ( 0) 
“ N w ( 0) 


V(t) 


1 

N w (r) 


^V(r w (T)) 


( 11 ) 


( 12 ) 


is the average of the potential energy over the random 
walker population, N w (t) is the instantaneous number 
of random walkers present in the population, and a is 
a proportionality constant, usually chosen to be 1/Ar. 
This (stabilization) procedure is imposed to keep the 
population size stationary in time throughout the course 
of simulation. Otherwise, without the stabilization, the 
population would explode or crash exponentially. Like¬ 
wise, the value of E le { fluctuates in time r in a stationary 
manner, and its time average can be used to estimate the 
ground state eigenenergy E 0 . Depending on the initial 
conditions for the random walker population, the sta- 
tionarity is achieved only at sufficiently long projection 
times. 

One should not be mislead by Eq. which gives the 
impression that the convergence of the present version 


of the DMC algorithm is defined by the maximum pro¬ 
jection time r max according to the rates with which all 
the unwanted contributions vanish exponentially (rela¬ 
tive to the contribution of the ground state \J/ 0 ) due to 
the factors e^ E °~ Ek ' >T ' aax . This would be the case if the 
ground state energy was estimated using the wavefunc¬ 
tion ^(r; r) explicitly. However, in the present algorithm 
it is estimated by averaging E le f, which fluctuates, and 
the convergence of the DMC energy estimate is rather 
defined by the statistical errors associated with this av¬ 
erage. 


Aside from the modifications associated with geometric 
constraints specific for the system in question, the present 
DMC algorithm exactly resembles the algorithm devel¬ 
oped by Anderson pji[2|, which was also implemented later 
by a number of groups (see, e.g. refs. HI 151 H51 [Ml 1521155)1 . 
Following ref. [14] we have augmented the branching step 
in order to prevent random walkers from escaping out 
of the basin of attraction associated with a chosen iso¬ 
mer. During the course of the simulation any walker that 
violates one of the imposed geometric constraints is elim¬ 
inated from the ensemble. Ideally, such constraints must 
be inexpensive to evaluate compared to the cost of the po¬ 
tential energy calculation, yet they should be effective in 
separating the isomer in question from the rest of the con¬ 
figuration space. The geometric constraints that we have 
imposed for simulations involving the water hexamer (see 
below) require calculations of the pairwise distances be¬ 
tween the oxygen atoms, distances of the oxygens from 
the center of mass, and the moments of inertia. Thresh¬ 
old values for the constraints are chosen as a compromise 
between restricting the configuration space only to the 
relevant basin of attraction and ensuring that all regions 
of the configuration space where the corresponding wave- 
function has an appreciable amplitude are accessible by 
the walkers. 


However, one must always keep in mind that there is 
only one true ground state, while the physical meaning 
of any other “ground state” of the system, subject to 
certain geometric constraints, must be established inde¬ 
pendently of the simulation algorithm. This problem is 
intimately related to the problem of defining an “isomer”, 
namely, a species that has a distinct structure, naturally 
associated with its own, local basin of attraction, a region 
in space that is unique enough to be separated from the 
other regions of configuration space by sufficiently large 
potential energy barriers. Alternatively, we can define 
an “isomer” as a state of the system, that once prepared, 
would exist for a sufficiently long time. Nevertheless, the 
actual stability of a quantum mechanical ground state 
restricted to a basin of attraction may not be trivial to 
establish, even approximately, without performing quan¬ 
tum dynamics simulations that are usually much more 
demanding than any ground state computation. 
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FIG. 1: Relative statistical errors for the ground state en¬ 
ergy estimates for H 2 O monomer and dimer as a function of 
1 / \JN w for a total projection time of 2.0 x 10 6 au. The er¬ 
rors were estimated by splitting the overall simulation into 10 
independent simulations. 
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A series of calculations with varying time step At and 
population size N w = N w (0) have been performed for 
each of the following: H 2 0, (H 2 0) 2 , D 2 0, and (D 2 0) 2 . 
Namely, for a fixed value of At = 10 au we ran six simu¬ 
lations with TV™ = 2xl0 2 , 4xl0 2 , 9xl0 2 , l.GxlO 3 , 3.6x 
10 3 , 1.96 x 10 4 , and for fixed N w = 1.96 x 10 4 , five sim¬ 
ulations with At = 2, 5, 10, 15, 20 au. Given At and 
N w , 10 independent DMC simulations were usually car¬ 
ried out. For each run a relatively short equilibration (if 
necessary) was performed before starting to accumulate 
the average energy (E le f). Consequently, the 10 energy 
averages were used to estimate the statistical error. The 
total projection time (i.e. including all 10 independent 
runs) for each At and N w was T max = 2 x 10 6 au. 

Clearly, the total numerical cost scales as 

numerical cost ~ T mSiX N w (13) 


and for given At and N w the statistical error should scale 
as ~ T ma x Because the total projection time is fixed, it 
follows from Eq. (13) that for a constant value of At, the 
numerical effort should grow linearly with population size 
N w . It is expected that for a fixed total projection time 
the statistical error should scale as ~ N w . This is 
demonstrated in Fig. [l] using the DMC results for the 
water monomer and dimer with T max = 2.0 x 10 6 au. 
Therefore, the statistical error should scale as 


statistical error ~ (N w T max ) 1 ^ 2 ~ (numerical cost) M 2 

Ml! 

That is, the statistical error is a function of the numerical 
cost only. Consequently, we assert that no computational 


time is gained through the use of small walker numbers, 
as for simulations with fewer random walkers, the projec¬ 
tion time must be increased accordingly to maintain the 
same statistical error. At the same time, the systematic 
bias in N w is reduced with increasing population size. 
This analysis suggests that in order to reduce the statis¬ 
tical error it is preferable to consider larger population 
sizes N w rather than longer simulation times T max . Of 
course, one should also take into account other consider¬ 
ations, such as limitations in the computer memory, and 
also the need to perform a sufficiently long equilibration, 
during which the averages are not accumulated. 

Fig. [2] and Fig. [3] show the dependence of the DMC 
energy estimates on the time step At, and the reciprocal 
walker number 1/N W , respectively. We observe that the 
time step bias for the ground state energy estimate at 
At = 10.0 au is noticeable. For example, for the water 
dimer it is ~ 0.03 kcal/mol. However, the estimate of the 
binding energy, 


Do ■= 2Ah 2 0 - -E(h 2 o) 2 (15) 

is hardly sensitive to the time step due to the nearly 
complete cancellation of the systematic errors (see the 
bottom left panel of Fig. 2|. At the same time, the 
population size bias for the binding energy estimate (the 
bottom left panel of Fig. [3]) is more evident and does 
not disappear completely. This behavior is attributed 
to the fact that the population of random walkers repre¬ 
sents the wavefunction and directly reflects the geometric 
composition of the system. For example, the population 
size bias is consistently stronger for the dimer than the 
monomer. In other words, the DMC energy estimate 
converges faster for the monomer than the dimer, which 
results in imperfect compensation of the systematic er¬ 
rors. Therefore, the failure of the biases for both systems 
to fully cancel appears in the binding energy curves as a 
residual bias (bottom left panel of Fig. [3]). Yet, for the 
isotope shift (the lowest right panel of Fig. [3]), 

SDo := Z?o(D 2 0) — D 0 (H 2 O) (16) 

the population size bias is significantly less pronounced. 
The systematic errors do cancel in this case, and the 
isotope shift converges to a value near SD 0 = 0.444 
kcal/mol. This weak dependence of the isotope shift on 
walker population indicates that the extent and behavior 
of the binding energy biases are similar for both types of 
isotopomers. 

Our results extrapolated to N w —> 00 are summarized 
in Table [I] Estimates of the statistical uncertainty in the 
DMC energies are at least an order of magnitude smaller 
than the systematic error (0.001 kcal/mol or lower) for 
the largest N w considered in Figs. [3] and [5] and, as such, 
are not included in the table. By extrapolating the ab¬ 
solute ground state energy estimates E 0 to the At — >• 0 
limit we could significantly improve their estimates. For 
example, for the water monomer H 2 0 we would obtain 
Eq = 13.18 kcal/mol, which coincides with the nurneri- 
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FIG. 2: DMC energies for the water dimer and monomer with N w = 1.96 x 10 4 walkers for five values of the time step At. 
Time step values range from 2.0 to 20.0 au. All calculations were run for a total projection time of 2.0 x 10 b au. Top Left: H 2 O 
dimer energies. Top Right: H 2 O monomer energies. Middle Left: D 2 O dimer energies. Middle Right: D 2 O monomer energies. 
Bottom Left: Binding energy Do for H 2 O and D 2 O dimers. Bottom Right: Isotope shift SDo for the dimer binding energies 
(note the scale). The “classical” energy (here and below) is computed by using the values of the corresponding potential energy 
minima. 
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FIG. 3: Same as Fig. [ 2 ] but as a function of the reciprocal walker number 1 /N w with fixed value of the time step At = 10.0 
au. Walker numbers range from 2.0 x 10 2 to 1.96 x 10 4 . 





















TABLE I: Classical ( D e ) and quantum (Do) estimates of the 
binding energies for the water clusters considered in this work. 
Eo stands for the ground state energy estimates. All the en¬ 
ergy estimates are obtained by extrapolating the DMC re¬ 
sults with Ar = 10.0 au to the N w —> oo limit. Experimen¬ 
tal binding energies for the dimer D^ pt ' were obtained from 
refs. Emm All energies are reported in kcal/mol. 
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ft ■ 
-‘-'mm 

Eo 

D e 

Do 

Dexpt. 

H 2 0 

0.00 

13.16 

- 

- 

- 

(H 2 0) 2 

-6.24 

21.80 

6.24 

4.53 

3.16 ±0.03 

(H 2 0) 6 -Cage 

-50.64 41.09 

50.64 37.89 

- 

(H 2 0)6-Prism 

-50.19 41.73 

50.19 

37.24 

- 

D 2 0 

0.00 

9.63 

- 

- 

- 

(d 2 o) 2 

-6.24 

14.29 

6.24 

4.97 

3.56 ± 0.03 

(D 2 0) 6 -Cage 

-50.64 

17.09 

50.64 40.61 

- 

(D 2 0)6-Prism 

-50.19 

17.69 

50.19 40.05 

- 


cally exact value obtained by diagonalizing the Hamilto¬ 
nian using a discrete variable representation. However, 
we emphasize that such extrapolation prior to taking the 
energy difference would be a bad idea, as in this case we 
would not be able to take advantage of the systematic 
error cancellations. Thus, while a bias does indeed exist 
for Eq, we determined that the time step error in Do is 
negligibly small for the dimer and is predicted to be so 
for the hexamer even with Ar = 10 au. 

We find that our DMC estimates of Do for both (^ 0)2 
and (D20)2 are about 1.5 kcal/mol higher than the ex¬ 
perimental values reported in refs. 13411351 which is a clear 
indication of the failure of the q-TIP4P/F PES in de¬ 
scribing accurately the energetics of small water clusters. 
Interestingly though, for the isotope shift in the binding 
energy the present result (~ 0.44 kcal/mol) happened to 
agree well with the experiment (~ 0.40 kcal/mol). 


Water Hexamer: Cage versus Prism. 

This section reports the results of DMC simulations for 
the water hexamer and its (D 2 0)6 isotopomer. Specif¬ 
ically, we have performed ground state calculations for 
two low-lying isomers: the cage and prism. The cage 
geometry corresponds to the global minimum of the q- 
TIP4P/F PES, while the prism local minimum is only 
0.45 kcal/mol higher. Originally, we also planned to con¬ 
sider the next-highest isomer in the sequence, the book, 
which lies about 1.2 kcal/mol higher than the cage mini¬ 
mum, but discovered that it is extremely unstable due to 
very low energy barrier separating it from the prism mini¬ 
mum. For example, in a classical Monte Carlo simulation 
a random walk starting in the book configuration would 
quickly jump into the prism minimum at temperatures as 
low as T ~ 5K, thus making the book isomer physically 
undetectable in a hypothetical experiment. Moreover, 
the quantum book isomer has even less physical mean¬ 
ing because including the quantum effects would only 
destabilize the system to a greater extent. In contrast 


to the book structure, the classical prism isomer appears 
to be stable up to temperatures as high as T ~ 60K, 
which makes the q-TIP4P/F prism geometry much less 
ambiguous to define and much more likely to be observed 
in a hypothetical experiment. However, the existence of 
relatively high potential energy barriers surrounding the 
prism isomer does not necessarily prevent the random 
walkers in the DMC simulation from leaking into other 
local minima, which happens often enough to make the 
problem nontrivial both conceptually and numerically. 
It is sufficient for only one random walker to escape into 
another energy minimum, where it may start replicating, 
thereby quickly causing a significant portion of the ran¬ 
dom walker population to represent the “wrong” geom¬ 
etry. As mentioned above, this problem can be circum¬ 
vented in principle by using importance sampling (see, 
e.g. refs. I23U26)) . However, these approaches may not be 
feasible to implement for systems as complex as water 
clusters, when an appropriate trial wavefunction must 
be very nontrivial to both determine and parametrize. 
Therefore, in the present context, we find the use of ge¬ 
ometric constraints less problematic. Intelligently cho¬ 
sen geometric constraints impose artificial barriers on the 
random walkers, preventing them from escaping out of 
the region defined by the particular isomer. Furthermore, 
the problem exists even in calculations of the true ground 
state, i.e., regardless of the relationship between the lat¬ 
ter and the isomer in question. In such a calculation the 
random walkers are initialized, usually but not always, in 
the global energy minimum (for the q-TIP4P/F PES, it 
is the cage minimum). Occasionally, one of the random 
walkers may still jump from the cage minimum to one of 
the local minima, where it has a chance to be replicated. 
Arguably, due to either physical or unphysical reasons, 
this migration of random walkers leads to their popula¬ 
tion being delocalized over energy minima representing 
different isomers. Consequently, for the present study, we 
find it necessary to impose geometric constraints even for 
DMC calculations involving the true ground state. 

The isomers of water clusters are usually identified by 
the arrangement of the oxygen atoms only; there might 
be several local minima, different only by the orientation 
of one or two water molecules, which are still identified by 
the same geometric motif (e.g. the prism). These min¬ 
ima are nearly isoenergetic, and they are often separated 
by low potential energy barriers. To this end, consider 
a configuration r = (iq,..., r^) (r, : e R 3 ) of K oxygen 
atoms (here K = 6) and a vector p p d[r] G 'R K ^ K ~ 1 ^ 2 
with elements defined by all the K{K — l)/2 pair dis¬ 
tances |r.; — r,j (i > j) arranged in descending order. 
The following constraint can then be deployed to pre¬ 
vent the random walkers from leaving the basin of at¬ 
traction defined by a reference configuration rh ef ) (e.g., 
corresponding to the prism minimum of the PES): 


Appd := 


y/2 


V K (K - 1) 


p pd [r] - p pd [r (ref) ] < hpd (17) 


where the full expression is the root-mean-square dis- 














9 


FIG. 4: The pair distance order parameter Ap p d (cf. Eq. (171) 
shown for several randomly chosen walkers as a function 
of projection time for two typical DMC simulations using 
N w = 8 x 10 5 with all the random walkers initialized in the 
prism minimum. Top: An unconstrained calculation in which 
the walkers started to leak out of the prism minimum into 
different regions of the PES at r ~ 10 5 au. Bottom: Same, 
except that geometric constraints were imposed to prevent 
walker leakage. 



Projection Time (au) 


placement of the instantaneous pair distances from those 
of the initial reference structure, and h . p d is an empir¬ 
ically chosen parameter. Note that the pair distance 
constraint 0 has two important properties: it is ro- 
tationally and translationally invariant and, secondly, its 
numerical evaluation is inexpensive, so it can be imple¬ 
mented at every MC iteration without a significant in¬ 
crease in the overall computational cost. However, due 
to its simplicity, it may not always be very effective. 

Consequently, we implemented two additional geomet¬ 
ric constraints, which follow from reasoning analogous to 
that used in formulating Eq. 0- The second constraint 
is based on the metric defined by the three principal mo¬ 
ments of inertia arranged in descending order within the 
vector p m i[r] € K 3 : 


A/9 mi := 


V3 


Pmi|r] - pn 


ef)l 


< K 


(18) 


Similarly, the vector p cm [r] £ R K with elements 
Yj — r cm | arranged in descending order contains the dis¬ 
tances between all K oxygen atoms and the center of 
mass r cm of the water cluster. Thus, the third constraint 
is defined by 


Ap cm 


1 

\fK 


r - pc 


ef)j 


< h c 


(19) 


The threshold values /i p a, h m ,, and h cm are established 
empirically by trial and error. For example, we discov¬ 
ered that the cage ground state calculation requires only 
the pair distances as a constraint with h pc j = 0.76 au, 


while for the prism ground state calculation all three con¬ 
straints are necessary with h p( j = 0.76 au, h m i = 125 
amu au 2 , and h cm = 0.66 au. 

Rg. a shows an example from two typical DMC simu¬ 
lations, without geometric constraints (top) and with all 
three constraints implemented using the prism as a ref¬ 
erence configuration. The quantity shown is that defined 
by Eq. 0 for several randomly chosen walkers as a func¬ 
tion of projection time r. The leaking out of the prism 
local minimum in the unconstrained simulation starts to 
occur at r ~ 10 5 au. 

Fig- a shows the DMC results for H 2 O hexamer prism 
and cage and their D2O isotopomers with fixed At = 10 
au and varying population size: N w = 4.0 x 10 4 , 6.0 x 10 4 , 
8.0 x 10 4 , 1.2 x 10 5 , 1.6 x 10 5 , 2.4 x 10 5 , 3.2 x 10 5 , and 
4.8 x 10 5 . The total projection time for most simulations 
was on the order of r max ~ 10 6 , which, as established ear¬ 
lier for the dimer, means that simulations with larger N w 
have smaller statistical errors. An examination of these 
plots reveals that for the population sizes considered, the 
bias in N w for the ground state energies of the hexamer is 
noticeable. For example, the change of the DMC energy 
estimate for the cage isomer is about 0.03 keal/mol when 
the population size is increased from N w = 3.2 x 10 5 to 
4.8 x 10 5 . Moreover, for N w = 4.8 x 10 5 , which is the 
largest walker number used in this work, the energy es¬ 
timate appears to be about 0.06 keal/mol higher than 
the value extrapolated to the N w —> 00 limit. At the 
same time, Fig. [5] shows that the biases in the absolute 
prism and cage energies follow a nearly identical trend re¬ 
gardless of which isotopomer is considered. Energies for 
(H 2 0) 6 appear to be slightly less converged than those 
for (D 20)6 at the same number of walkers and projec¬ 
tion times (i.e., statistical errors are smaller for the “less 
quantum” (D 2 0)6). The slopes of the prism and cage 
energies as a function of 1 /N w are virtually equivalent 
(even up to small features), which indicates that the bias 
is not dramatically affected by variations in the spatial 
arrangement of the monomers. As a consequence, upon 
taking the energy difference between the prism and cage 
isomers, the bias is substantially reduced resulting only 
in a weak dependence on N w . Visual inspection of the 
energy curves enables us to conclude that our DMC es¬ 
timates of the ground state energy differences, 

A E = Eo (prism) — Eq( cage) (20) 

are accurate to about 0.01 keal/mol, when sufficiently 
large walker numbers (i.e., N w > 2.0 x 10 5 ) are consid¬ 
ered. Additionally, the small oscillations in 1/N W that 
still remain in the DMC energy estimates correlate be¬ 
tween (H 2 0) 6 and (D 2 0) 6 , which makes us believe that 
these oscillations are real. 

Considering now the isotope shift (as shown in the 
lower left panel of Fig. [ 5 ]), 

^isotope ■— AFh 2 0 AEb 2 o (21) 

we conclude that further cancellations of the systematic 
errors occur for this quantity, such that the result seems 
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1/N 1/N 

w w 


FIG. 5: DMC energies for the water hexamer with At = 10 au as a function of the reciprocal walker number 1 /N w . The 
walker numbers used in the simulations range from N w = 4.0 x 10 4 to 4.8 x 10°. Top Left: Ground state energies Eo for 
(H20)6 cage and prism. Top Right: Same, for (D20)6 cage and prism. Bottom Left: Prism-cage energy differences for 
(H20)6 and (D20)6 isotopomers A E = .Eo(prism) — _E 0 (cage). Bottom Right: Isotope shift for prism-cage energy differences 
<5isoto P e = A£h 2 o — AB D2 o (note the scale). 


to reach a plateau near isotope = 0.055 ± 0.01 kcal/mol 
at N w = 1.2 x 10 5 . However, compared to this small 
value of the isotope shift, the bias for walker numbers 
N w < 1.2 x 10 5 is noticeable, and at N w ~ 5.0 x 10 4 
isotope even changes sign. 

Conclusions 

We have provided an assessment of the DMC method 
of Anderson. In particular, we have analyzed the bias 
with respect to the time step At and the walker number 
N w through application to the water monomer, dimer, 
and hexamer along with their D 2 O isotopomers. We have 
determined that the use of relatively small walker num¬ 
bers N w is unprofitable because longer projection times 
An ax are then required to maintain the same statistical 
error. Consequently, the best-case scenario is to choose 


larger values of N w , as both the statistical and systematic 
errors decrease with increasing N w . 

The calculations undertaken in this work for the hex¬ 
amer used walker numbers up to N w = 4.8 x 10 5 with 
total projection times on the order of r max ~ 10 6 , which 
for At = 10 au required ~ 10 5 x N w potential en¬ 
ergy evaluations per isomer. Nevertheless, the absolute 
ground state energies did not level off to a well-defined 
plateau. As demonstrated here, obtaining accurate esti¬ 
mates for the binding energies does require high accuracy 
for the absolute values of the ground state energies. In¬ 
deed, the population size bias for the water dimer does 
not cancel with that of the water monomer, for which 
the bias is negligible. However, this is not an issue for 
the time step bias, as the systematic errors for the dimer 
and monomer are removed almost completely upon de¬ 
termination of the binding energy. At the same time, the 
biases undergo cancellation for interesting properties of 
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the water hexamer such as the prism-cage energy differ¬ 
ence and isotope shift. Therefore, it is possible to obtain 
accurate values for these quantities in practice. How¬ 
ever, without carrying out a systematic study, as done 
in this work, one cannot assume that the errors would 
fortuitously cancel. Also, the q-TIP4P/F PES renders 
our calculations less expensive by orders of magnitude 
compared to presumably more accurate ab initio -based 
surfaces EHHsII] or even true ab initio surfaces. Should we 
use one of these surfaces, this publication in 2015 would 
hardly be possible. 

Interestingly, the quantum effect for the cage-prisnr 
energy difference, <5q M := AE H2 o - A£ c i assica i = 0.20 
kcal/mol, is relatively small, while the isotope shift 
isotope = 0.055 kcal/mol is even smaller. Both these 
results are consistent with an earlier prediction in which 
the Self Consistent Phonons method was applied to the 
water hexamer using several different water mo dels [36]. 
In light of this conclusion, we cannot ignore the DMC 
result of ref. [TS] regarding the water hexamer (albeit us¬ 
ing the WHBB PES of ref. |29|) , which indicated that the 
cage and prism bound state energies are nearly the same, 
E 0 (prism) « E 0 (cage), while the classical cage-prism en¬ 
ergy difference for this PES is A£’ c i ass i ca i « 0.4 kcal/mol. 


The maximum number of random walkers N w = 1.6 x 10 s 
and the total projection time r max = 8.0 x 10 5 au used 
in ref. ns give rise to statistical errors noticeably larger 
than those in the present study. However, the most im¬ 
portant difference between the present DMC study and 
that of ref. M is that no geometric constraints were im¬ 
plemented in the latter, while we found the geometric 
constraints to be of crucial importance in preventing the 
random walker population from spilling out of the poten¬ 
tial energy minimum corresponding to a particular iso¬ 
mer. If such a procedure is not implemented, one should 
keep his or her fingers crossed in the hope that during 
the course of the DMC simulation the walker population 
does not spread over several local minima (corresponding 
to the “wrong” isomers), thereby resulting in an incorrect 
energy estimate. 
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